Method and apparatus for determining viscoelastic parameters in tissue

ABSTRACT

Described herein are a method and apparatus for determining viscoelastic parameters of a tissue. A vibration signal is applied to the tissue and displacements at a plurality of locations within the tissue are measured at a plurality of times. The viscoelastic parameters of the tissue, including elasticity and viscosity, can then be determined by fitting a finite element model of the tissue to the vibration signal and the measured displacements and by solving for the viscoelastic parameters of the model. A value for density of each element of the model is selected and the absolute values for the viscoelastic parameters of each of the elements in the model is determined. Alternatively, the difference in relaxation-times between two locations within the tissue can be determined from the difference in phases of the strains at the two locations.

CROSS REFERENCE TO RELATED APPLICATION

This application claims priority under 35 U.S.C. §120 to a U.S.non-provisional application having application Ser. No. 12/611,736,filed Nov. 3, 2009, which claims priority under 35 U.S.C. §119(e) to aU.S. provisional application having Ser. No. 61/193,185, filed Nov. 3,2008, the entireties of which are herein incorporated by reference.

TECHNICAL FIELD

The present invention relates to a method and apparatus for determiningone or more viscoelastic parameters in tissue.

BACKGROUND

Viscoelastic parameters of tissue include the elasticity and viscosityof the tissue. These viscoelastic parameters can vary between healthyand diseased regions of tissue. For example, cancerous breast tumors areapproximately eight times stiffer than benign lesions while benignlesions are approximately four times stiffer than normal breast tissue(J. Ophir, I. Cespedes, B. Garra, H. Ponnekanti, Y. Huang and N. Maklad,“Elastography: ultrasonic imaging of tissue strain and elastic modulusin vivo,” European Journal of Ultrasound, vol. 3, no. 1, pp. 49-70,January 1996). Therefore, a doctor examining a patient for a tumor willtypically press down on various regions of the patient's body, searchingfor relatively stiff (highly elastic) regions that could signify thepathology.

Viscosity is also an important mechanical property when assessing tissuehealth. For example, cancerous lesions with a high blood vesselconcentration have a higher viscosity than surrounding, healthy tissue(R. Sinkus, M. Tanter, T. Xydeas, S. Catheline, J. Bercoff, and M. Fink,“Viscoelastic shear properties of in vivo breast lesions measured by MRelastography,” Magnetic Resonance Imaging, vol. 23, no. 2, pp. 159-165,2005).

Tissue can be modeled as being composed of a series of interconnectedelements, with each element being composed of a mass (m) coupled to aparallel arrangement of a spring having a spring coefficient (k) and adamper having a damping coefficient (b). The spring is used to model theelasticity of the tissue, which is reflected in the tissue's resilienceto an applied external force. The damper is used to model the viscosityof the tissue, which is reflected in the tissue's internal friction andrelaxation effects. The mass is used to model the density of the tissue.When relatively small applied forces are applied and relatively smalldisplacements result, the force-velocity relationship of the dampers islinear, as is the force-displacement relationship of the springs. Allobjects with mass exhibit a linear relationship between force andacceleration. One, two or three dimensional arrays of elements are alsopossible. For an example of a two dimensional array, see FIG. 10 of U.S.patent application having Ser. No. 10/963,795 and publication number2005/0119568 to Salcudean et al., the entirety of which is hereinincorporated by reference.

To estimate relative values for elasticity and viscosity in tissue,vibrations can be applied to the tissue and the resulting tissuedisplacement observed. This can be done over a one, two or threedimensional region and using a variety of imaging modalities; ultrasoundand magnetic resonance imaging (MRI) are both modalities that have beenused. See, for example, U.S. Ser. No. 10/963,795 to Salcudean et al.,and Emre Turgay, Septimu Salcudean, and Robert Rohling, “Identifying themechanical properties of tissue by ultrasound strain imaging,”Ultrasound in Med. & Biol., vol. 32, no. 2, pp. 221-235, 2006, theentirety of which is herein incorporated by reference. By “relativevalue”, it is meant a value for elasticity or viscosity at a certainlocation within the tissue that is expressed relative to anotherlocation within the tissue. A “relative value” is in contrast to an“absolute value”, which is a value for elasticity or viscosity that isexpressed independently from any other location within the tissue. Forexample, an absolute value for viscosity is 10 Pas, whereas a relativevalue for viscosity is 10·x Pas, where “x” is a baseline value forviscosity measured elsewhere in the tissue.

However, prior art methods are often slow and cannot image a tissueregion in real-time, can require expensive machinery, or can onlydetermine relative values of viscoelastic parameters.

Consequently, there exists a need for a method and system that addressesone or more deficiencies in the prior art.

BRIEF DESCRIPTION OF THE DRAWINGS

In the accompanying drawings, which illustrate one or more exemplaryembodiments:

FIG. 1 depicts a one-dimensional mass-spring-damper model of tissue;

FIG. 2 depicts an ultrasound probe applying a vibration signal to thetissue;

FIG. 3 depicts the ultrasound probe applying the vibration signal to thetissue and illustrates the axial, lateral, and elevational planes of thetissue;

FIG. 4 is a schematic view of a hand-held actuator having the ultrasoundprobe affixed thereto;

FIG. 5 is a block diagram of a system for determining viscoelasticparameters in tissue, according to one embodiment;

FIG. 6 is a flowchart describing how viscoelastic parameters in the formof elasticity, relaxation-time, and viscosity can be determined;

FIG. 7 is a flowchart describing how relaxation-times can be determinedaccording to another embodiment;

FIG. 8 is a graph of relaxation-time for the tissue vs. frequency,determined according to the embodiment of FIG. 7;

FIG. 9 is a flowchart describing how viscoelastic parameters can beiteratively solved, according to another embodiment;

FIG. 10 is a flowchart depicting how viscoelastic parameters can benon-iteratively solved, according to another embodiment;

FIG. 11 depicts a two-dimensional region of tissue divided into varioussegments;

FIGS. 12A, 12B, 12C, and 12D depict graphs of transfer functionsdetermined according to the embodiment of FIG. 7; and

FIGS. 13A, 13B, 14A, 14B, 15A, and 15B depict images of viscosity andelasticity determined in accordance with the aforementioned embodiments.

DETAILED DESCRIPTION

Typical viscoelastic parameters of a region under study that aremeasured or desirable include elasticity (E), viscosity (η), andrelaxation-time (T), which is defined as T=η/E. For the purposes of theembodiments discussed herein, the region under study is human tissue,such as breast tissue. The embodiments described herein are directed ata method and apparatus for determining one or more of elasticity,viscosity and relaxation-time. Two approaches are employed to accomplishthis.

Referring now to FIG. 1, the tissue can be modeled as interconnectedelements each having a mass (m) 12 coupled to a parallel arrangement ofa spring 14 having a spring coefficient (k) and a damper 16 having adamping coefficient (b). The first approach involves assuming thatvibration of the tissue is sufficiently slow to be able to ignore theinertial effect of the mass 12, thereby allowing the mass 12 to beremoved from the model. Removing the mass 12 from the model allows thetissue region to be modelled as a series of interconnected spring-damperelements. Ignoring the mass 12 of each element is reasonable for avibration signal of low frequency (e.g.: less than 50 Hz). A “vibrationsignal” is any externally applied stimulus to tissue that causes tissuedeformation that can be measured by means of an imaging modality, suchas ultrasound. Due to the vibration signal, tissue is deformed as afunction of time.

The tissue model is depicted in FIG. 1. In FIG. 1, the vibration signalis applied to the tissue, which results in an external force f(t) beingapplied to the tissue. The tissue is divided into elements, with eachelement extending between two nodes. From the measured axialdisplacements at the two nodes surrounding any particular element, thestrain experienced by that element can be calculated. Followingcalculation of strain at any two elements, the transfer functionrelating the strains experienced by those two elements can becalculated. For example, the transfer function H₂ ¹(ω) relates thestrains between elements 1 and 2 such that ε₁(ω)=H₂ ¹(ω)·ε₂(ω). Inaddition to relating strains between any two elements, a transferfunction can also be calculated that relates a strain experienced by oneelement to a stress experienced by a second element, or between adisplacement experienced by a first element to a displacementexperienced by a second element. This first approach is hereinafterreferred to as the “Transfer Function Approach”.

The second approach involves assuming that tissue deformation can bemodeled by a linear viscoelastic finite element model. Instead ofignoring the mass as is done in the Transfer Function Approach, thedensity is presumed to be close to 1,000 kg/m³, which is the density ofwater, since water is the main constituent of most soft tissues. If moreaccurate and local a priori information of the density in differentparts within the tissue are available (e.g.: from textbooks), thatinformation can be incorporated into the model as well. Such a model canbe linear or non-linear in the remaining unknown parameters. The linearequations can be solved iteratively or non-iteratively; the equationsthat are not linear can be solved iteratively. This second approach ishereinafter referred to as the “Finite Element Model Approach” or “FEMApproach”.

Both the Transfer Function Approach and the FEM Approach are implementedusing the same hardware.

FIG. 2 depicts a handheld actuator 24 that is used to apply a vibrationsignal to tissue 28. An operator holds the actuator 24 in his or herhand 22. The actuator 24 is connected to an ultrasound imaging system(not shown in FIG. 2) through a cable 18 which drives a motor inside theactuator 24 and which produces motion and, consequently, the vibrationsignal. The vibration signal can be a filtered wide-band signal coveringmultiple frequencies (i.e.: be composed of multiple, superimposedfrequency components). Alternatively, the vibration signal can be asinusoidal waveform composed of only a single frequency component. Themoving part of the actuator 24 is connected to an ultrasound probe 26which is in direct contact with or otherwise coupled to the tissue 28,such as by a water-based gel and produces mechanical waves in the tissue28. The boundary conditions 30 of the tissue 28 may be known, unknown orpartially known. While the vibration signal is applied to the tissue 28via the ultrasound probe 26, the ultrasound probe 26 receivesradio-frequency (RF) signals and transmits the signals to the ultrasoundmachine through a cable 20 for processing and display.

The actuator 24 is shown in more detail in FIG. 4 and can be attached tothe ultrasound probe 26 to produce a mechanical stimulus in the tissue.A fastener 13 tightens the ultrasound probe 26 to a vibrating stage 2.The ultrasound probe 26 is connected to the ultrasound machine by thecable 20. A static stage 1 is held by the operator and is connected tothe motor 3 by another fastener 4. The motor 3 is connected to anoblique extension 8 by a motor shaft 7. An actuation system transformsmotor rotational motion into reciprocating translational motion througha variable angle swash plate 9, in an arrangement similar to that of aswash plate pump. The swash plate angle can be adjusted manually with ascrew 10 while the mechanical stimulus is applied. The swash plate 9 isconnected to the vibrating stage 2 by a hinge 11. A linear guide 15ensures linear motion of the vibrating stage 2 with respect to thestatic stage 1.

A block diagram of an apparatus for determining viscoelastic parametersin tissue is depicted in FIG. 5. The hand-held actuator 24 andultrasound probe 26 are integrated into a larger system of hardwareconfigured to receive ultrasound data from the tissue 28, determineviscoelastic parameters from the ultrasound data, and the display theviscoelastic parameters to the operator. The probe 26 is pressed overthe tissue 28, which is limited by physical boundaries defining theboundary conditions. Included in the system of FIG. 5 is an ultrasoundmodule 500, which is part of the ultrasound imaging system, thattransmits ultrasound signals to and receives ultrasound data from theultrasound probe 26 and controls the ultrasound pulse formation. Theultrasound module 500 transmits data to a parameter estimation module502 that performs processing steps including, but not limited to, motionestimation, the estimation of relaxation-time and power-law parametersas used in the Transfer Function Approach, the estimation of viscosityand elasticity using finite element modeling as used in the FEMApproach, and post-processing of parameters such as viscosity andelasticity. Post-processing of parameters can include feature extractionand classification of datasets. The transformed data generated in theparameter estimation module is transmitted to a display module 506,which generates image maps of the desired estimated parameters andrepresents these on a display, which can include a computer monitor. Anactuation module 504 drives the mechanical excitation in the hand-heldactuator 24. The actuation module 504 transmits a mechanical excitationsignal through a controller and an amplifier to drive the actuator 24.The actuation module 504 performs tasks including, but not limited to,waveform processing and data pipelining of the mechanical excitationsignal, performance monitoring of the mechanical excitation, and agraphical user interface (GUI) that allows the user to adjust parametersincluding, but not limited to, the frequency (for the case of a singlefrequency excitation), bandwidth or frequency range (for a the case ofexcitation at multiple frequencies), amplitude and waveform shape of themechanical excitation signal. The parameter estimation module 502, theultrasound module 500, the actuation module 504 and the display module506 can form different modules running on hardware such as theultrasound machine, a general purpose computer, or any suitable dataprocessor coupled to a memory having stored thereon one or more of themodules 500, 502, 504 and 506. Alternatively, any one or more of themodules 500, 502, 504 and 506 can be stored on a computer readablemedium such as random access memory and any form of disk based media.Also depicted in FIG. 5 are an amplifier and controller communicativelycoupled to the actuator 24, which are part of a driver circuit thatdrives the actuator 24 to generate the vibration signal.

If ultrasound is used as the imaging modality for viscoelasticitymeasurement of the tissue 28, the ultrasound probe 26 is placed on thesurface of the tissue 28 and acquires radio-frequency (RF) data alongseveral A-lines (straight lines extending from the surface of the probe26) in one, two or three dimensions (1D, 2D or 3D). Therefore, eachpoint in the tissue will have a displacement in 3D which can beprojected into axial (along an A-line), lateral (perpendicular to axialin the imaging plane) and elevational (perpendicular to the imagingplane) directions; these directions are depicted in FIG. 3. Axial,lateral and elevational strains are the spatial derivative of suchdisplacements. By applying a Fourier transform or a similar frequencytransform to the time-domain displacement or strains, the displacementor strain phasors can be calculated at a desired frequency (ω). Forexample, the phasor of the axial displacement of a certain location atfrequency w can be obtained by calculating the Fourier transform of thatdisplacement signal at ω.

Referring now to FIG. 6, there is depicted a flowchart describing howviscoelastic parameters in the form of elasticity, relaxation-time, andviscosity can be determined according to the Transfer Function Approach(referred to as “Estimation method 1” in FIG. 6) and a 2D version of theFEM Approach (referred to as “Estimation method 2” in FIG. 6). At block600, ultrasound RF signals are acquired from the tissue 28. If theTransfer Function Approach is to be employed, axial motion along onedimension of the tissue 28 is derived from the ultrasound RF signals atblock 602 and the Transfer Function Approach is employed at block 604,resulting in elasticity and relaxation-time. If a 2D version of the FEMApproach is to be employed, one or both of axial motion and lateralmotion is derived from the ultrasound RF signals at blocks 602 and 606,and the FEM Approach is employed at block 608, resulting in absolutevalues of elasticity and viscosity for each of the elements of thefinite element model. Notably, a 3D version of the FEM Approach can alsobe used, in which axial, lateral and elevational motion are considered.

The Transfer Function Approach and the FEM Approach are now described inturn.

The Transfer Function Approach

When mass 12 is assumed to be zero, the tissue 28 can be modeled as aplurality of elements connected in series, with each element composed ofa parallel arrangement of a spring and a damper; this is called a Voigtmodel, with each element being called a Voigt element. The relationshipbetween stress (σ) and strain (ε) of a Voigt element is σ(t)=Eε(t)+η{dotover (ε)}(t), where the dot is the time derivative operator, E isYoung's modulus (also known as elasticity), and η is the viscosity. Inthe frequency domain, stress and strain are related by {circumflex over(σ)}=(E+jωη){circumflex over (ε)}, where the hat operator represents afrequency transform, or phasor. Therefore, by defining relaxation-timeas τ=η/E, the phase difference between the stress and strain is∠{circumflex over (σ)}−∠{circumflex over (ε)}=tan⁻¹(ωτ).

Given that relaxation-time is on the order of milliseconds, atrelatively low frequencies (e.g.: less than 50 Hz for typical softtissues), the phase difference can be approximated as ∠{circumflex over(σ)}−∠{circumflex over (ε)}=ωτ. While strain can be measured at all theVoigt elements throughout the tissue, stress values cannot be measuredat Voigt elements at interior points of the tissue. Furthermore,accurate measurement of forces at the area of contact between theultrasound probe 26 and the surface of the tissue requires 2D arrays offorce sensors, which may not be practical. However, by assuming a onedimensional spring-damper tissue model, the stress in the tissue can beassumed to be equal everywhere along any given one dimensional line.Consequently, ωτ₁=ωτ₀+(∠{circumflex over (ε)}₀−{circumflex over (ε)}₁)where terms having an index of 0 correspond to a first location on theline and terms having an index of 1 correspond to a second location onthe line. The difference in relaxation-times between two locations onthe line can therefore be determined from the difference of the phase ofthe strains between those two locations.

Referring now to FIG. 7, there is depicted a flow chart illustrating howrelaxation-times of the tissue are calculated in a region of tissuemodeled using N Voigt elements. At block 700, axial displacements alongone A-line that are caused by application of the vibration force aremeasured. From the measured axial displacements, the axial strain iscalculated at the various elements at block 702. At block 704, thetime-dependent strain of one of the Voigt elements (ε_(r)) is used as abaseline value to obtain relative relaxation-time measurements. If thevibration signal includes only a single frequency component (block 706),then at block 716 the phase difference between the strains at thereference element ε_(r) and the strain at any other element ε_(i) can becalculated as ∠{circumflex over (ε)}_(r)−∠{circumflex over (ε)}_(i).From this, the relative relaxation-time of ε₁ relative to ε_(r) isτ_(i)=−(∠{circumflex over (ε)}_(i)−∠{circumflex over (ε)}_(r))/ω, ascalculated in block 718.

If the vibration signal applied to the tissue 28 is composed of a singlefrequency component, the phasors can be calculated at that frequency. Ifthe vibration signal has multiple frequency components, the phasors canbe calculated at all of the plurality of frequencies. The axialdisplacement or axial strain phasors can be used with a 1D model oftissue in order to measure the relaxation-time. Transfer functionsbetween the displacement or strain signals at two locations can becalculated at a desired frequency by computing the ratio between theirFourier transforms at that frequency. By “transfer function”, it ismeant any transformation that yields the ratio of the phasors at the twolocations; a transfer function can be determined by calculating theFourier transform of the two signals and dividing them at a desiredfrequency. Alternatively, a transfer function can be determined bydividing the cross-spectral density of the two signals by theauto-spectral density of the reference signal.

The transfer functions can be calculated at a single frequency to beable to measure the relaxation-time at that frequency. Alternatively,the transfer functions can be calculated at a plurality of frequenciesto measure the relaxation-time at each of those frequencies or averagedto obtain an average of the relaxation-times at those frequencies. Ifthe vibration signal includes multiple, superimposed frequencycomponents, then instead of proceeding to block 716 from block 706,transfer functions H_(r) ^(i)(ω) are calculated at block 708. H_(r)^(i)(ω) is the transfer function between strains measured betweenelements r and i such that ε_(i)(ω)=H_(r) ^(i)(ω)·ε_(r)(ω).Consequently, the difference in relaxation-time between elements i and ris

${{{\tau_{i}(\omega)} - {\tau_{r}(\omega)}} = {\frac{- 1}{\omega}\angle \; {H_{r}^{i\;}(\omega)}}},$

as executed at block 710. If the relaxation-time at element r is known,the absolute relaxation-time of element i can be calculated at block712. In an alternative embodiment (not depicted) in which the stresses(σ) at elements i and r are known and can be assumed equal, the absoluterelaxation-time at any given element, having a strain ε, can becalculated as

${\tau = {\frac{- 1}{\omega}\angle \; {H_{\sigma}^{ɛ\;}(\omega)}}},$

where H_(σ) ^(ε)(ω) is the transfer function between the stress and thestrain signals.

If the relaxation-time is known at one location along a line undergoingdynamic displacement in a specimen through independent calibration orrheometry tests, the relaxation-time can be determined at allfrequencies for all the other points on that line. Once such profilesare obtained for the relaxation-times, given the non-Newtonian behaviorof soft tissues, a power-law model can be fit to the relaxation-timecurves. At block 714, power law parameters (T and α) of every element inthe tissue are calculated. The power-law model is used to interpolatethe mechanical characteristics of the material at a desired frequency orto use the power-law model parameters for further analysis, featureextraction or classification. Non-Newtonian fluids are those fluids forwhich the viscosity changes with frequency. Generally, in fluids andsoft, viscoelastic solids, the viscosity decreases by a power-law as thefrequency increases. FIG. 8 shows the power-law decrease of therelaxation-time for a few tissue mimicking materials, including a softbovine-skin gelatin phantom (12 wt. % in water)—Gelatin A; a hardgelatin phantom (18 wt. % in water)—Gelatin B; and a wet PVA sponge atthree states: completely soaked in water—Sponge C; Partiallysqueezed—Sponge B; and severely squeezed—Sponge A. The slope ofrelaxation-time vs. frequency is an additional mechanical characteristicof a viscoelastic material that physicians can use to furthercharacterize different types of tissue.

In alternative embodiments (not depicted), other viscoelastic parameterscan be derived from the transfer function H_(r) ^(i)(ω) or H_(σ)^(ε)E(ω). Above, transfer function phase is used to measurerelaxation-time; additionally, Young's modulus can also be obtained fromreciprocal of the magnitude of the transfer function between the strainand the stress (|H_(σ) ^(ε)(ω)|⁻), or relative estimation can beobtained from the reciprocal of the magnitude of the transfer functionbetween the two strains at two locations (|H_(r) ^(i)(ω)|⁻). In oneembodiment, the viscosity can be estimated by multiplying an estimate ofYoung's modulus and the relaxation-time at every pixel location. Inother embodiments, a power-law model can be fitted to the viscosity as afunction of frequency, with the goal of using the frequency-dependencyof the viscosity estimate as another power-law parameter for tissuedifferentiation.

In all of the aforedescribed embodiments, the relaxation-time parameter(τ) or the power-law parameters (T and α) of every block in the tissuecan be displayed on the screen to form an image, where the pixel valuesof the image represent viscoelastic parameters and are calculated inreal-time or offline. Physicians can use these images to identify normaland diseased tissues or different types of tissue.

In various embodiments, displacements are used to determine transferfunctions. In other embodiments, measurements of strain or displacementsin other directions may be used, such as the lateral and elevationaldirections. For example, the transfer functions can be computed betweentwo different axial displacements, two different lateral displacements,two different lateral strains or between any two points of othercomponents of the displacements or strains.

Parameters derived from the phase of the transfer function, such asrelaxation-time and viscosity, can be displayed as an image where pixelvalues represent the value of the parameter. In some embodiments, animage representing the viscoelastic parameters can be superimposed on asecond image such as a conventional ultrasound, MRI, or CT image, toprovide context for interpreting the viscoelastic parameters.

The FEM Approach

Introducing mass or density into finite element model (FEM) formulationsof viscoelastic materials enables the estimation of the absolute valuesof the elasticity and viscosity parameters rather than the relativevalues only even when stress or force is unknown. Such quantitativemeasurement can be achieved even without having any measurement of theapplied force. An FEM model that takes into account the viscous andinertial effects is a discrete representation of the wave equation in aviscoelastic medium. Therefore, quantitative measurement of theviscoelastic parameters can be performed based on knowledge orassumptions about the density of the elements. If the density of all ofthe elements in the model are known based on a priori knowledge such asfrom a textbook, that knowledge can be used in the model to estimate theabsolute value of the viscoelastic parameters. In most cases, thisdensity distribution is not known exactly, but it is reasonable toassume that water density dominates in soft tissue. Thus the density canbe assumed to be 1000 kg/m³. Furthermore, in the absence of forcemeasurements, quantitative measurement of other parameters such as wavespeed, kinematic viscosity and relaxation-time can be achieved. Withaccurate and quantitative characterization of viscosity at differentfrequencies in this manner, a power-law trend can also be decipheredfrom the viscosity and/or relaxation-time at different frequencies andthe power-law parameters can be estimated and used to identify differenttissue types.

In order to measure the elasticity and viscosity parameters by fittingthe measured displacement data to a finite element model, time-domaindisplacements are measured in 1D, 2D or 3D from a region of interest intissue. The displacement phasors are calculated from the time domaindisplacement measurements. The region of interest is meshed by severalinterconnected elements. Adjacent elements are connected to each otherat their vertices, called nodes. The displacements and forces may onlybe measured on or applied to the nodes of the mesh. FIG. 11 depicts a 2Dfinite element model mesh 1100 having elements 1102 and nodes 1104.

Each element has an elasticity, a viscosity and a density property. Thedensity is assumed to be a certain value as described above. Thedisplacement phasors are measured at the nodes 1104. Subsequently, thecomplex-valued displacement phasors are fitted to a frequency-domainrepresentation of the dynamic finite element model and quantitativeviscoelastic parameters of all the elements are measured by solving aninverse problem.

The FEM model can represent a 2D or 3D region. In the 2D case, axial andlateral displacements or a subset thereof are required to solve the FEMinverse problem. In the 3D case, axial, lateral and elevationaldisplacements or a subset thereof are required to solve for theviscoelastic parameters of the model. Unlike the Transfer FunctionApproach, the FEM Approach accounts for 2D and 3D displacementmeasurements. Furthermore, in the FEM Approach, the viscoelasticparameters can be measured quantitatively without needing to know theapplied forces, whereas absolute measurements are not possible withoutknowing applied forces or stresses in the Transfer Function Approach.

The axial and lateral displacement phasors or time-domain displacementdata can be used with a 2D finite element model to measure viscosity andelasticity in a 2D region of interest in the tissue 28. The axial,lateral and elevational displacement phasors or time-domain displacementdata can be used with a 3D finite element model to more accuratelymeasure viscosity and elasticity than in a 2D region of interest. Ifvolumetric measurement of displacements are available, the viscoelasticparameters can be similarly measured for a volume, using a 3D FEM model.

In the FEM Approach, the tissue 28 is divided into elements 1102 forwhich viscoelastic parameters are derived based on a fit to a model thatsatisfies input and output signals. The input signals are the nodalboundary conditions and the output signals are the displacements of theinterior nodes 1104 of the model. The nodal boundary conditions are thedisplacements at the boundaries of the tissue 28 or the external forcesto the model. These external forces are usually not measured, althoughit is possible in specific cases to measure the pattern of forcesapplied by the actuator 24. To represent the viscoelasticcharacteristics of the tissue 28 the following transient FEM model isused, with the displacement vector u(t) describing the motion at all thenodes 1104 in response to applied external nodal forces f(t):

Ku(t)+B{dot over (u)}(t)+Mü(t)=f(t)  (1)

where K is the stiffness matrix which depends on the elements'elasticity values, concatenated into a vector E; B is the viscousdamping matrix which depends on the elements' viscosities, concatenatedinto a vector η; and M is the mass matrix which depends on the elements'densities, concatenated into a vector p. Assuming lumped masses at thenodes 1104 of the FEM model results in M being a diagonal matrix.

For a linear FEM model, each entry of the K matrix, called the stiffnessmatrix, is a linear combination of the elasticities of the elements 1102in the model.

Similarly, the entries of B, called the damping matrix, and M, calledthe mass matrix, are linearly dependent on the element viscosities anddensities, respectively. Any technique as is known to a person skilledin the art can be used to construct the K, B and M matrices in an FEMmodel.

In one embodiment, the damping matrix can be constructed for everyelement 1102 from the stress-strain relationship in a linearviscoelastic medium. In this embodiment, the damping matrix of everyelement 1102 is linearly dependent on the viscosity of the tissue 28 atthat element 1102. The global damping matrix B is then obtained bycombining the individual damping matrices for each element 1102. In thismanner, the entries of the damping matrix B will be linearly dependenton the viscosity of the elements 1102.

To derive B according to this embodiment, note that the stress tensor(σ_(ij)) in the Cartesian coordinate system for a viscoelastic materialcan be written as follows,

σ_(ij)=λε_(kk)δ_(ij)+2με_(ij)+λ′{dot over (ε)}_(kk)δ_(ij)+2μ′{dot over(ε)}_(ij),  (1a)

The strain tensor is denoted by and its time derivative by {dot over(ε)}_(ij). Tensor notation is used; thus ε_(kk)=ε₁₁+ε₂₂+ε₃₃ and theKronecker delta δ_(ij)=1 if i=j and 0 if i≠j. The subscript indicescorrespond to orthogonal directions in the Cartesian coordinate system.λ and μ are the Lamé constants and λ′ and μ′ are the viscositycharacteristic parameters. In most fluids and rubber-like materials, thebulk viscosity which is defined as κ=λ′+2μ′/3 is approximately zero,thus λ′=−2μ′/3. The vector notation of (1a) is:

σ=Cε+C′{dot over (ε)},  (1b)

where vectors σ=[σ₁₁ σ₂₂ σ₃₃ σ₂₃ σ₃₁ σ₁₂]^(T) and ε=[ε₁₁ ε₂₂ ε₃₃ ε₂₃ ε₃₁ε₁₂]^(T) contain the six components of the isotropic stress and straintensors in the general 3D setting. In (1b), C is the material elasticitycharacteristic matrix which can be derived from (1a). C′ is the materialviscosity characteristic matrix which can be derived from (1a) with thezero bulk viscosity assumption:

$\begin{matrix}{C^{\prime} = {\frac{2}{3}{{\mu^{\prime}\begin{pmatrix}2 & {- 1} & {- 1} & 0 & 0 & 0 \\{- 1} & 2 & {- 1} & 0 & 0 & 0 \\{- 1} & {- 1} & 2 & 0 & 0 & 0 \\0 & 0 & 0 & 3 & 0 & 0 \\0 & 0 & 0 & 0 & 3 & 0 \\0 & 0 & 0 & 0 & 0 & 3\end{pmatrix}}.}}} & \left( {1c} \right)\end{matrix}$

While the above notation is for a 3D FEM problem, for a 2D case thematrices C and C′ can be reduced by assuming plane-stress orplane-strain as known to the skilled in the art. For each element 1102,the viscous damping matrix B can be produced by integrating (1b) overthe entire element with appropriate shape functions.

Given that the entries of K are also linearly dependent on theelasticity of the elements 1102, the FEM equations can be established asa system of equations that are linear in the viscoelastic parameters.Beneficially, when constructing the damping matrix according to thisembodiment, a linear system of equations can be obtained which can besolved, as described in greater detail below.

Conventional methods of constructing the damping matrix B do not yield alinear system of equations. For example, in one method of constructingthe damping matrix known as the Rayleigh damping, the damping matrix isa linear combination of the stiffness and mass matrices, B=ξ₁K+ξ₂M,where ξ₁ and ξ₂ are constant coefficients that describe the dependenceof B on K and M. In this representation, the viscosity of the elementsare not incorporated into the damping matrix and the FEM system ofequations will not be linear in the viscosity parameters of the elements1102. Therefore, the preferred approach is not to use the standarddamping matrix formulation, which leads to a non-linear dependence onparameters, but rather employ the linear stress-strain relationship in aviscoelastic medium to construct the damping matrix.

By taking the Fourier transform of Equation (1), the linear FEMformulation of the tissue in the frequency domain is

(K+jωB−ω ² M)û={circumflex over (f)}  (2)

where hatted variables denote Fourier transforms and j is the imaginaryunit. Noting that û and {circumflex over (f)} are complex-valuedvectors, equation (2) can be decomposed into its real and imaginarycomponents:

(K+jωB−ω ² M)(û ^(r) +jû ^(i))={circumflex over (f)} ^(r) +j{circumflexover (f)} ^(i)  (3)

where superscripts r and i represent the real and imaginary components,respectively, and j is the imaginary unit.

By decomposing û and {circumflex over (f)} into their real and imaginaryparts, equations (2) and (3) can be written as:

$\begin{matrix}{{A\begin{bmatrix}{\hat{u}}^{r} \\{\hat{u}}^{i}\end{bmatrix}} = \begin{bmatrix}{\hat{f}}^{r} \\{\hat{f}}^{i}\end{bmatrix}} & (4)\end{matrix}$

where the superscripts r and i denote the real and imaginary componentsof a vector, respectively. A is defined as follows:

$\begin{matrix}{A = {\begin{bmatrix}{K - {\omega^{2}M}} & {{- \omega}\; B} \\{\omega \; B} & {K - {\omega^{2}M}}\end{bmatrix}.}} & (5)\end{matrix}$

If the viscoelastic parameters of the FEM model are known, the problemof solving for the displacement vector u as a result of a known forcevector f and certain displacement boundary conditions is called theforward problem. If the displacements are known and Equations (4) and(5) are solved for the elasticity and viscosity of the elements, it isthe inverse problem that is solved. Note that the elasticity andviscosity of the elements have linear relationships with the entries ofthe matrices K and B, respectively.

Two approaches can be employed to solve the inverse problem. In a firstapproach, Equations (4) and (5) can be solved in an iterative fashion,without assuming linearity in the unknown viscoelastic parameters. In asecond approach, the linearity of Equations (4) and (5) in the unknownviscoelastic parameters can be exploited by solving Equations (4) and(5) with either an iterative or non-iterative method. Note that Equation(4) describes a system of equations that is linear in the viscoelasticparameters. In the first approach, the linearity in parameters is notexploited in finding a solution. In the second approach, the linearityin parameters is exploited by using standard methods to solve a linearsystem of equations.

Solving Equations (3) and (4) Iteratively, without Assuming Linearity

To be able to solve the forward problem to obtain the displacements atthe interior nodes 1104, given a distribution of the viscoelasticparameters, the displacement boundary conditions can be incorporatedinto Equation (2), as is known to persons skilled in the art. Thisprocedure involves changing the externally applied forces on the righthand side of Equation (2), as well as the corresponding rows in the K, Band M matrices on the left-hand side of Equation (2), to account for thedisplacement at the boundary, instead of the force. This way, a newleft-hand side matrix and right hand side vector are obtained,subscripted by c, as in the following equation:

$\begin{matrix}{{G_{c}\begin{bmatrix}{\hat{u}}^{r} \\{\hat{u}}^{i}\end{bmatrix}} = \begin{bmatrix}{\hat{f}}_{c}^{r} \\{\hat{f}}_{c}^{i}\end{bmatrix}} & (6)\end{matrix}$

The subscript c indicates that the corresponding matrix or vector hasbeen manipulated to replace the equations corresponding to the unknownexternal forces by the corresponding displacement boundary conditions.G_(c) is defined as follows:

$\begin{matrix}{G_{c} = \begin{bmatrix}{K_{c} - {\omega^{2}M_{c}}} & {{- \omega}\; B_{c}} \\{\omega \; B_{c}} & {K_{c} - {\omega^{2}M_{c}}}\end{bmatrix}} & (7)\end{matrix}$

Without displacement boundary conditions, solving Equation (4) for thetissue elasticity and viscosity in the region of interest is anill-posed problem. After applying the displacement boundary conditions,Equation (6) is well-posed and can be solved to obtain the displacementsin the interior region.

As an example, we assume that there is only one non-zero entry (f₀) inthe force vector which has to be replaced by the corresponding knowndisplacement boundary condition (u₀=U0^(r)). Therefore, assuming thatthe first column of G in (4) is [g₁₁, g₂₁ . . . g_(ng1)]^(T), where ngis the number of rows in G, Equation (4) will change to:

$\begin{matrix}{{\begin{bmatrix}1 & 0 & \ldots & 0 \\0 & \; & \; & \; \\\vdots & \; & G_{s\; 11} & \; \\0 & \; & \; & \;\end{bmatrix}\begin{bmatrix}{\hat{u}}^{r} \\{\hat{u}}^{i}\end{bmatrix}} = \begin{bmatrix}{U\; 0^{r}} \\{{- g_{21}}U\; 0^{r}} \\\vdots \\{{- g_{{ng}\; 1}}U\; 0^{r}}\end{bmatrix}} & \left( {7\; e} \right)\end{matrix}$

where G_(s11) is a matrix of G obtained by removing the first row andthe first column of G. The resulting left-hand side matrix is G_(c) in(6) that is computed for this example and the resulting right-hand sidevector is the

$\quad\begin{bmatrix}{\hat{f}}_{c}^{r} \\{\hat{f}}_{c}^{i}\end{bmatrix}$

vector that is calculated for this example.

In solving the inverse problem, the optimal values of Young's moduli andviscosities of the elements 1104 are determined iteratively to fit themodel according to the observed displacement vectors. To solve theinverse problem iteratively and determine optimal values for Young'smodulus and viscosity, a cost function can be used to measure deviationof predicted displacements with estimated parameters from measureddisplacements. The cost function can be defined according to Equation(8).

$\begin{matrix}{{\Phi (p)} = {{\frac{1}{2}{{{{\hat{u}}^{r}(p)} - {\hat{u}}_{0}^{r}}}^{2}} + {\frac{\alpha}{2}{{{{\hat{u}}^{i}(p)} - {\hat{u}}_{0}^{i}}}^{2}}}} & (8)\end{matrix}$

In this equation, p=[E^(T) η^(T)]^(T) is a 2m parameter vector composedof elasticity and viscosity parameters, considering that there are melements in the FEM model. The parameter α is a weight determining therelative importance of the imaginary and real components of thedisplacement in the inverse problem. α can simply be equal to one. TheEuclidean norm of a vector is denoted by ∥.∥.

A number of cost functions similar to (8) can be utilized to solve formodel parameters, and the function in (8) has been used only as anexample of a particular embodiment.

The parameter vector p can be updated using an iterative procedureincluding, but not limited to, the modified Gauss-Newton method, theLevenberg-Marquardt algorithm or any other method that finds a set ofparameters that minimizes a cost function such as Equation (8). In thepresent embodiment, the parameter vector is updated iteratively as,p_(k+1)=p_(k)+Δp_(k), k=1, 2, . . . , where k is the iteration index.The update vector Δp_(k) is to be calculated from:

H _(k) Δp _(k) =−J _(k) ^(T) Δû _(k),  (9)

where J_(k) is the Jacobian, H_(k) is the Hessian and Δû_(k) is definedas follows:

$\begin{matrix}{{{\Delta \; {\hat{u}}_{k}} = {\begin{bmatrix}{\Delta \; {\hat{u}}_{k}^{r}} \\{\Delta \; {\hat{u}}_{k}^{i}}\end{bmatrix} = \begin{bmatrix}{{\hat{u}}_{k}^{r} - {\hat{u}}_{0}^{r}} \\{{\hat{u}}_{k}^{i} - {\hat{u}}_{0}^{i}}\end{bmatrix}}},} & (10)\end{matrix}$

where, for simplicity û_(k) denotes û(p_(k)). By applying the chain ruleto calculate the derivative of (6), the following results:

$\begin{matrix}{{{{\frac{\partial G_{c}}{\partial p_{j}}\begin{bmatrix}{\hat{u}}^{r} \\{\hat{u}}^{i}\end{bmatrix}} + {G_{c}\left\{ J_{j} \right\}}} = \begin{bmatrix}\frac{\partial{\hat{f}}^{r}}{\partial p_{j}} \\\frac{\partial{\hat{f}}^{i}}{\partial p_{j}}\end{bmatrix}},} & (11)\end{matrix}$

where {J_(j)} denotes the j^(th) column of the Jacobian matrix. Notethat p_(j) can be either elasticity or viscosity of an element,therefore ∂G_(c)/∂p_(j) can be computed from (7) as follows:

$\begin{matrix}{\frac{\partial G_{c}}{\partial p_{j}} = \left\{ \begin{matrix}{\begin{bmatrix}\frac{\partial K}{\partial e_{j}} & 0 \\0 & \frac{\partial K}{\partial e_{j}}\end{bmatrix},} & {p_{j} = e_{j}} \\{{\omega \begin{bmatrix}0 & {- \frac{\partial B}{\partial\eta_{j}}} \\\frac{\partial B}{\partial\eta_{j}} & 0\end{bmatrix}},} & {p_{j} = \eta_{j}}\end{matrix} \right.} & (12)\end{matrix}$

Given that G_(e) is known, Equation (11) can be solved for {J_(j)} andthus the Jacobian matrix can be constructed.

To solve for Δp_(k) in equation (9), the Hessian at iteration k iscalculated using (13).

H _(k) =J _(k) ^(T) J _(k)+λ_(k) I _(2m×2m),  (10)

In (10), I_(2m×2m) is an identity matrix with m being the number ofelements, and λ_(k) is a small regularization factor that makes theH_(k) positive definite.

In alternative embodiments, a spatial filter can be applied to theelasticity and viscosity distributions at every iteration to make thesolution of the problem smooth and less sensitive to the displacementnoise. As a result, the elasticity or viscosity parameter of eachelement 1102 will be a weighted sum of its own value and the values ofthat parameters in the adjacent elements 1102. A linear filter can thusbe constructed in the form of a sparse matrix that contains the requiredweights for each element 1102. The filter can be convolved with theupdated distribution of the parameter at every iteration to conduct theoptimization toward a smooth solution.

Referring now to FIG. 9, there is depicted a method for solvingEquations (6) and (7) iteratively, which can be executed on anultrasound machine. Beginning at block 900, a user measuresdisplacements phasors û₀ at a given frequency ω. At block 902, aninitial distribution for p is assumed and, subsequently, at block 904,û(p) is determined based on the value of p assumed at block 4902. Theinitial guess for p can be found by setting all elasticity parametersand all the viscosity parameters of the different elements in the FEMmodel equal to each other. For example, all the elasticity values can beinitialized to 20 kPa and all the viscosity values can be initialized to20 Pas. At block 908 the processor can check whether lateraldisplacements were measured in the interior of the region of interest ofthe tissue 28 or not. If lateral displacements were measured, theprocessor proceeds to block 912, where p is updated using anoptimization method such as that described in Equation (9). However, ifthe lateral displacements were not measured, the columns of the Jacobianthat correspond to those lateral displacements can be removed and areduced Jacobian can be obtained (block 410). This reduced Jacobian canthen be used to calculate the Hessian as described above, and to solvefor Δp_(k) in Equation (9) using only the axial portion of Δû_(k), and pcan be updated in block 912. If lateral displacements are not measured,φ(p) will be measured based on the differences in the axialdisplacements only. At block 914, φ(p) and its gradient are calculated.The value of φ(p) is an indication of how close to the optimal solutionthe current parameter value is, while its gradient indicates how slowlyor quickly the procedure is converging. Whether or not one of theseconvergence criteria are met is checked at block 916. If the convergencecriterion is not met, the procedure is repeated beginning at block 904.If the convergence criterion is met, then the method finishes at block918.

Solving Equations (3) and (4) Assuming Linearity in the Parameters,Either Non-Iteratively or Iteratively Deriving Ap=b

When assuming linearity, a series of equations are obtained from thefinite element model in the form of Ap=b, where p represents theviscoelastic properties, i.e. the elasticities and viscosities of allthe elements in the model. The description below provides the specificprocedure used in this embodiment to obtain a linear system of equationsin the form of Ap=b from the finite element model in Equation (3).

Starting from Equation (3), the entries of the stiffness matrix K, areformed as a linear combination of the elasticity parameters of theelements:

$\begin{matrix}{{K_{l,q} = {\sum\limits_{z = 1}^{m}\; {c_{l,q}^{z}E^{z}}}},} & (14)\end{matrix}$

where c_(l,q) ^(z) are constant coefficients that are determined by thetype of mesh and geometry. Furthermore, m is the number of elements andK_(l,q) refers to an entry on row l and column q of the stiffnessmatrix. Similarly, the global damping matrix and global mass matrix arelinearly dependent on the viscosity and density parameters,respectively.

Given that the entries of K are linearly dependent on the elasticity ofthe elements, the entries of vectors Kû^(r) and Kû^(l) from (3) are alsolinearly dependent on the elasticity of the elements. For instance, thel^(th) entry of Kû^(r) is as follows:

$\begin{matrix}{{\left\lbrack {K{\hat{u}}^{r}} \right\rbrack^{l} = {\sum\limits_{q = 1}^{dn}\; {K_{l,q}{\hat{u}}_{q}^{r}}}},} & (15)\end{matrix}$

where n is the number of nodes in the model and d is the number ofdimensions in the problem which is equal to 2 for 2D and 3 for 3D cases.Substituting the expression for K_(l,q) from (14) into (15) andreformulating that as a vector-vector product:

$\begin{matrix}\begin{matrix}{\left\lbrack {K{\hat{u}}^{r}} \right\rbrack^{l} = {\sum\limits_{q = 1}^{dn}\; {\sum\limits_{z = 1}^{m}\; {c_{l,q}^{z}E^{z}{\hat{u}}_{q}^{r}}}}} \\{{= {C_{l}^{{Er}^{T}}E}},}\end{matrix} & (16)\end{matrix}$

where E is the elasticity vector and C_(l) ^(Er) is a vector defined asfollows:

$\begin{matrix}{{C_{l}^{{Er}^{T}} = {\sum\limits_{q = 1}^{dn}\; {{\hat{u}}_{q}^{r}c_{({l,q})}^{T}}}},} & (17)\end{matrix}$

Given that one element of the vector Kû^(r) can be calculated through avector-vector multiplication, the entire Kû^(r) vector can be computedby performing a matrix-vector multiplication as follows:

Kû ^(r) =C ^(Er) E,  (18)

where the superscripts E and r indicate that this coefficient matrixcorresponds to the elasticity parameters and was constructed using thereal part of the displacement vector. Similarly, for all of thecombinations in (3):

Kû ^(r) =C ^(Er) E,

Kû ^(i) =C ^(Ei) E,

Bû ^(r) =C ^(ηr)η,

Bû ^(i) =C ^(ηi)η,

Mû ^(r) =C ^(ρr)ρ,

Mû ^(i) =C ^(ρi)ρ,  (19)

where η and ρ are the vectors containing the viscosities and densitiesof the elements. By substituting (19) into (3) and splitting the realand imaginary parts into separate equations, the following results:

$\begin{matrix}{{\begin{bmatrix}C^{Er} & {{- \omega}\; C^{\eta \; i}} \\C^{Ei} & {\omega \; C^{\eta \; r}}\end{bmatrix}\begin{bmatrix}E \\\eta\end{bmatrix}} = {{{\omega^{2}\begin{bmatrix}C^{\rho \; r} \\C^{\rho \; i}\end{bmatrix}}\rho} + {\begin{bmatrix}{\hat{f}}^{r} \\{\hat{f}}^{i}\end{bmatrix}.}}} & (20)\end{matrix}$

The rightmost term in (20) is the external force vector of the FEM. Notethat the entries of the force vector are zero everywhere except for afew nodes 1104 to which the external excitation is applied. Theequations that correspond to non-zero forces can be removed from therows of the matrices in (20). The following model is thus obtained whichrepresents a system of equations that is linear in the parameters p:

Ap=b  (21)

where,

$\begin{matrix}{{A = \begin{bmatrix}C^{Er} & {{- \omega}\; C^{\eta \; i}} \\C^{Ei} & {\omega \; C^{\eta \; r}}\end{bmatrix}_{R}},} & (22) \\{{p = \begin{bmatrix}E \\\eta\end{bmatrix}},} & (23) \\{{b = {{\omega^{2}\begin{bmatrix}C^{\rho \; r} \\C^{\rho \; i}\end{bmatrix}}_{R}\rho}},} & (24)\end{matrix}$

where the subscript R denotes that the corresponding matrices arereduced in size by eliminating the rows that correspond to theexternally applied forces.

As an example, assume that there is only one non-zero entry (f₀) in theforce vector which has to be eliminated from the equations. Assumingthat f₀ is the first entry in the force vector, on the right-hand sideof (20), the first rows of the matrices (20) and consequently the firstequation in the system of equations described by (20) should be removed.The remaining system of equations can be described as in (21).

Equation (21) can be solved for parameters p, if the unknown externalforces are removed from the equations. Therefore, the rows in A and bthat correspond to unknown external forces are removed. To calculate Aand b based on row elimination for unknown forces, two cases should beconsidered. In the first case, if a node 1104 has a zero external forcecomponent in a specific direction, the corresponding row in the matricesin (20) should not be changed or removed. In the second case, if a node1104 has a non-zero external force component at a specific direction,the corresponding row should be removed from equation (20) to be able tosolve (21).

In solving (21), A is a known matrix which depends on the mesh, regionof interest, frequency and measured displacements. The variable p is theunknown parameter vector comprising elasticity and viscosity parametersof all of the elements. The vector b is a known vector which depends onthe mesh, region of interest, frequency, measured displacements anddensity of the elements 1102. As discussed above, the density can beassumed to be a known value, and approximately equal to that of water(1000 kg/m³) for all elements. Alternatively, the density is not assumedto be constant across all elements but is assumed to be known based on apriori knowledge, such as from textbooks on the typical properties of atissue type.

In an alternative embodiment, the unknown force entries can be insertedinto the unknown parameter vector, p, instead of being removed. In thiscase, one column should be added to the A matrix for each unknown forceentry. Therefore the problem can be solved simultaneously for theabsolute values of the viscoelastic parameters as well as the unknownforces, assuming that there are more equations than unknowns in thesystem. For example, assume that there is only one unknown force entry,f₀, which is real-valued and applied to the first node. In this caseequation (21) holds, while equations (22-24) will be changed to:

$\begin{matrix}{{A = \begin{bmatrix}C^{Er} & {{- \omega}\; C^{\eta \; i}} & {- 1} \\C^{Ei} & {\omega \; C^{\eta \; r}} & 0\end{bmatrix}},} & \left( {22\; e} \right) \\{{p = \begin{bmatrix}E \\\eta \\f_{0}\end{bmatrix}},} & \left( {23\; e} \right) \\{{b = {{\omega^{2}\begin{bmatrix}C^{\rho \; r} \\C^{\rho \; i}\end{bmatrix}}\rho}},} & \left( {24\; e} \right)\end{matrix}$

where, −1 in (22e) is a scalar while 0 is a vector of zeros. By using(22e) and (24e) in (21), the absolute values of the viscoelasticparameters and unknown force in (23e) can be measured.

In another embodiment, the finite element equations of the viscoelasticmodel are analyzed in the time domain (1) to measure the quantitativevalues of the parameters. In this embodiment, a time sequence ofdisplacements is measured and used in the time-domain representation ofthe dynamic finite element model. Therefore, an equation similar to (19)can be derived with time-domain displacement matrices instead of real orimaginary displacement phasors. As a result, a similar equation to (20)but with larger matrices can be obtained and the same procedure can beused to solve the resulting system of equations. Absolute values of theelasticity and viscosity parameters of all of the elements can beobtained by solving an inverse problem.

It is possible to solve for p in (21) if there exist more equationscompared to the number of unknown parameters. Therefore, afterelimination of the unknown external forces, the number of unknownparameters must be less than the number of equations in (21).

Considering a 2D FEM mesh with m×n nodes, comprising rectangularelements, the number of all the available equations is equal to the twotimes the number of nodes (2mn). The number of unknown viscoelasticparameters for all the elements is two times the number of elements(2(m−1)(n−1)=2mn−2(m+n−1)). Therefore, there are 2(m+n−1) more equationsthan unknowns in the case where there is one set of measurements at onefrequency. In this case, in order to be able to solve the inverseproblem, the number of unknown entries in the force vector cannot bemore than 2(m+n−1). A similar calculation can be performed for the 3DFEM model and for different types of elements.

The number of equations in the system of equations can be increased byincorporating more displacement measurements. In one embodiment, moredisplacement measurements at two or more frequencies increase the numberof equations. For example, by having measurements at two frequenciesinstead of only one frequency, and incorporating them into the A matrixin (21), the number of rows in A will double. As a result, more robustestimates of the parameters can be obtained. Incorporating multiple setsof displacements at multiple frequencies into A is also useful to reducethe sensitivity to the measurement noise in the displacements.Therefore, even if there are sufficient number of equations to solve theproblem, displacements at multiple frequencies can be used to constructA and obtain a more robust estimate of the viscoelastic parameters.

If the elimination of the unknown entries from the force vector resultsin insufficient number of equations to solve for all of the parameters,other element reduction techniques can be performed as described below.

Parameters derived from the finite element model such as viscosity andelasticity can be displayed as an image where pixel values represent thevalues of the particular parameter of interest. In some embodiments, animage representing the viscoelastic parameters can be superimposed on asecond image to give additional information, and to give a context forinterpreting the viscoelastic parameters. The second image can be, forexample, an ultrasound image, a magnetic resonance image, or a computedtomography image.

Solving Ap=b

Equation (21) can be solved in different ways, as described below.

(a) Solving the Unconstrained Problem

If there are no additional equality or inequality constraints on thevariables, Equation (21) can be solved by employing several algorithmsto solve a minimum least squares problem. In one embodiment of thisapproach, A^(T) can be pre-multiplied by the two sides of (11) to obtainthe following solution for p:

p=(A ^(T) A)⁻¹ A ^(T) b  (25)

Here, p can be computed using techniques known to persons skilled in theart.

In another embodiment of this approach, (21) can be solved usingiterative strategies such as the conjugate gradient method.Alternatively, (21) can be solved using singular value decomposition,Cholesky factorization or QR factorization.

In another approach, instrumental variables can be used to solve (21).In one embodiment of this approach, an instrument matrix A_(v) can beused to pre-multiply the two sides of (21). This gives:

A _(v) ^(T) Ap=A _(v) ^(T) b  (26)

Here again, p can be computed as is known to persons skilled in the art.

The instrument matrix should have low or zero correlation with themeasurement noise embedded in A. The instrumental variables techniquecan be applied in a recursive way to improve the estimation of p withtime. The instrument matrix is generated in the same way as A, but usinga different set of measurements. In one embodiment, the instrumentmatrices can be computed from the measurements at a different timeinstance. Alternatively, the instrument matrices can be computed fromthe measurements at a different frequency.

The instrumental variable technique has been applied in other fieldssuch as process control identification, but has not been used before inthe estimation of mechanical or viscoelastic parameters of tissue.Advantageously, it provides for the accumulation of measurements inorder to determine more reliable parameter estimates. Without its use,noise in the measurements can cause the estimates to severelyunder-estimate the actual parameter values.

(b) Solving the Constrained Problem Linear Programming (Non-Iterative)

In another approach, a linear programming technique can be used to solve(21).

In one embodiment of this approach, a convex optimization problem can beset up in which the minimization functional is zero or any constantscalar. In such a linear programming problem, (21) can be applied as aconstraint to the problem. Other linear equality and inequalityconstraints can be added to the optimization problem while maintainingits “linear program” formulation. For example, the elasticity andviscosity parameters of the model can be constrained within a range thatis known for soft tissue. Furthermore, a limit can be imposed on thespatial variation of the parameters. For example, an inequalityconstraint can be assumed on a linear combination of the parameters

Quadratic Programming (Iterative)

Another approach to solve (21) is to set up the problem as a quadraticprogram and solve it using one of the various methods to solve aquadratic program. The quadratic program can be defined as follows:

$\begin{matrix}{{\min\limits_{p}{\frac{1}{2}p^{T}{Hp}}} + {z^{T}p}} & (27)\end{matrix}$

Linear equality and inequality constraints on the parameters of themodel, as described above for the case of a linear program, can be addedto the optimization problem while maintaining its “quadratic program”formation. In one embodiment of this approach, a cost function can bedefined as ∥Ap−b∥², and can be expanded as:

∥Ap−b∥ ² =p ^(T) A ^(T) Ap−2b ^(T) Ap+b ^(T) b  (28)

Knowing that b^(T)b is a scalar and does not affect the minimization,from (21) and (28), it can be deduced that:

H=2A ^(T) A,  (29)

z=−2A ^(T) b,  (30)

Therefore, the following quadratic program is obtained:

$\begin{matrix}{{{\min\limits_{p}{p^{T}A^{T}{Ap}}} - {2\; b^{T}{Ap}}},{{subject}\mspace{14mu} {to}\mspace{14mu} \left\{ \begin{matrix}{{{A_{eq}p} = b_{eq}},} \\{{lb} \leq {A_{i}p} \leq {ub}}\end{matrix} \right.}} & (32)\end{matrix}$

where A_(eq) and A_(i) are matrices that comprise the coefficients ofthe equality and inequality constraints, respectively. b_(eq), lb and ubare vectors defining the equality constraints, lower and upper bounds,respectively.

In another embodiment of this approach, the instrumental variable methodcan be used to set up the minimization problem as a quadraticprogramming method with or without additional constraints on theparameters.

Elimination of the Parameters in the Finite Element Model Using EqualityConstraints

A limited number of elements 1102 can be assumed to have the sameviscoelastic parameters. This way, the number of unknown parameters inthe problem will be reduced and the inverse problem can be solved forthe viscoelastic parameters. Assuming that certain elements 1102 havethe same parameters is known as imposing “equality constraints”.

In one embodiment, a limited number of adjacent elements 1102 can beassumed to have the same viscoelastic parameters. This condition can beapplied to the aforementioned constrained optimization methods as a fewequality constrains to equate the viscoelastic parameters to thoseadjacent elements. In this way, the number of unknown parameters in theproblem can be reduced and the inverse problem can be solved more easilywithout significantly affecting the results for parameters removed fromthese adjacent elements.

In another embodiment of this approach, the tissue 28 is divided intoone or more sets of segments, such that the viscoelastic parameters forall of the elements 1102 inside each segment are constant. In FIG. 11,for example, shaded area 1106 represent a segment in which all elements1102 have the same viscosity and elasticity values. The number ofunknown parameters will only be twice the number of segments, as all theelements 1102 in the segment share common elasticity and viscosityparameters. This way, the number of the unknown parameters will bereduced and solving the FEM inverse problem will be more efficient.

Alternatively, this can be implemented by assuming several equalityconstraints between elements. As an example of this approach, if a tumoris detected or suspected by a clinician, its boundaries can bedelineated and a region of interest can be defined to enclose the tumor.The entire region of interest can be divided into two segments withdifferent viscoelastic properties. Thus the absolute elasticity andviscosity can be measured inside the tumor and in the background tissue.In this manner, a characterization of the tumor properties is possiblein an efficient way and the processing unit can measure and display theparameters at a faster rate, enabling real-time imaging.

According to another embodiment, homogeneity of the parameters in theentire region can be assumed, thus measuring the absolute values of onlytwo parameters, such as elasticity and viscosity, of the homogeneousregion, without knowing the forces.

Referring now to FIG. 10, there is depicted a flowchart explaining howthe parameter estimation module 502 can solve the linear system ofequations Ap=b. At block 1000, the operator acquires data by measuringdisplacements within the tissue 28 using the ultrasound probe 26.Following data acquisition, at block 1002 the parameter estimationmodule 502 generates an indexed notation of the finite element modelproblem, as described in Equations (20) and (21).

Subsequently, at block 1004, the parameter estimation module 502constructs a system of linear equations in the form Ap=b as described inEquations (21) through (24). The density of the elements of the finiteelement model is assumed to be a certain value, for example the densityof water, at block 1006. The system of linear equations is subsequentlysolved at block 1008, optionally taking into account any constraints atblock 1010. At block 1012, the parameter estimation module 502 hassolved the equations.

Application to Other Imaging Modalities

Aside from ultrasound, other imaging modalities, such as magneticresonance imaging (MRI) and computed tomography (CT), can also be usedto acquire axial and lateral displacements of tissue in response to amechanical excitation.

EXAMPLES Example 1 Using the Transfer Function Approach to DeriveRelaxation-Time: Numerical Results

A medium with four Voigt elements has been simulated with the elasticityand viscosity as shown in FIGS. 12A and 12B, respectively. A wide-bandmotion has been applied to the first node and the time-dependentdisplacements and strains have been calculated. The magnitude and phaseof the transfer functions between the stress and strains are shown inFIGS. 12C and 12D. Note that the low-frequency behavior of the magnitudeof the spectrum correlates with the elasticity of the medium, while theslope of the phase depends on the local relaxation-times.

Example 2 Using the FEM Approach to Determine Viscoelastic Properties ina PVC Phantom

The FEM Approach was used to reconstruct elasticity and viscosity of acubic PVC phantom (6 cm×5 cm×3 cm) with a hard circular inclusion. Theimages were produced as follows: an ultrasound linear probe acquired RFdata at 98 frames per second from the top of the phantom forapproximately 1 second while an actuator was applying 20 Hz vibration tothe phantom from below, having a peak-to-peak amplitude of less than 1mm. Radiofrequency data were collected from a 4 cm×4 cm region in 2D.Two dimensional displacements were calculated and their real andimaginary parts were computed at 20 Hz and the quadratic programmingmethod was used to estimate the absolute values of the elasticity andviscosity using a finite element model. FIGS. 13A and 13B show thereconstructed elasticity and viscosity maps of the PVC phantom with ahard circular inclusion. The absolute value of the elasticity is closeto the palpable elasticity of the phantom.

Example 3 Using the FEM Approach to Determine Viscoelastic Properties ina Gelatin Phantom

A tissue mimicking phantom (38 mm×40 mm×25 mm) has been constructed with12% (by weight) bovine skin gelatin in water. A small piece of apolyvinyl alcohol (PVA) sponge was embedded in the phantom. Based onindependent rheometric measurements, the background gelatin had aYoung's modulus and a viscosity of approximately 15 kPa and 15 Pasrespectively at 10.7 Hz. The sponge was approximately three times harderand 15 times more viscous than the background at that frequency.

The phantom was placed on an actuator and ultrasound RF data werecaptured using a linear probe positioned at the top of the phantom. Awide-band displacement excitation (1-30 Hz) with a Gaussian amplitudedistribution with standard deviation of 44 μm was applied to the bottomof the phantom. The internal motion of the phantom was estimated using across-correlation based algorithm. Frequency analysis was performed onthe displacements and the real and imaginary parts at 10.7 Hz wereselected to solve the inverse problem by minimizing a convex functionaliteratively. The estimated contrast of the elasticity and viscositybetween the background and inclusion are in agreement with therheometric measurements of the viscoelastic properties as noted above.

The elasticity and viscosity distributions measured in the phantom aredepicted in FIGS. 14A and 14B, respectively.

Example 4 Using the Transfer Function Approach to DetermineRelaxation-Time in a Phantom

A DC motor was used to apply a wideband excitation (3-30 Hz) to thephantom from the bottom, with less than 0.3% overall strain. RF data wasacquired at 98 frames per second over a 40 mm imaging window from thetop of the phantom. The phantom was constructed from 12% bovine skingelatin. A hard gelatin inclusion (18% by weight) was placed on the leftside of the phantom while a small PVA sponge inclusion was embedded onthe other side of the phantom. Axial displacements were calculated andthe transfer functions between the strains at every line were computed.The asymptotic magnitude of the transfer function between 3 Hz and 10 Hzwas used to generate the relative elasticity image and the phase of thetransfer functions at 15 Hz was used to generate the relativerelaxation-time image of the phantom. The relative values between thebackgrounds and inclusion are in agreement with independent rheometricmeasurements. The relative elasticity and relaxation-time images aredepicted in FIGS. 15A and 15B, respectively. The images indicate thatwhile both inclusions are harder than the background material, only thesponge inclusion has a significantly higher relaxation-time compared tothe surrounding gelatin.

While particular embodiments have been described in the foregoing, it isto be understood that other embodiments are possible and are intended tobe included herein. It will be clear to any person skilled in the artthat modifications of and adjustments to the above embodiments, notshown, are possible.

1. A method for determining viscoelastic parameters of a tissue, themethod comprising: (a) applying a vibration signal to the tissue; (b)measuring displacements at a plurality of locations within the tissue inresponse to the vibration signal at a plurality of times; and (c)determining the viscoelastic parameters of the tissue by fitting afinite element model of the tissue to the vibration signal and themeasured displacements and solving for viscoelastic parameters of themodel, wherein a value for density of each element of the model isselected and absolute values for one or more of the viscoelasticparameters for each of the elements of the model are solved.
 2. A methodas claimed in claim 1 wherein the measured displacements are expressedin a frequency domain as phasors having a real and an imaginarycomponent.
 3. A method as claimed in claim 1 wherein the viscoelasticparameters that are solved include both elasticity and viscosity.
 4. Amethod as claimed in claim 1 wherein the value for density of eachelement is selected to be the density of water.
 5. A method as claimedin claim 3 wherein the viscoelastic parameters further comprise one ormore of relaxation-time and wave speed.
 6. A method as claimed in claim1 wherein the vibration signal consists of one frequency component.
 7. Amethod as claimed in claim 1 wherein the solving for viscoelasticparameters comprises generating a system of equations comprising aviscosity and elasticity of each of the elements of the model andwherein the system of equations is treated as being linear in theviscosity and elasticity of each of the elements.
 8. A method as claimedin claim 7 wherein the viscoelastic parameters are unconstrained.
 9. Amethod as claimed in claim 7 wherein the viscoelastic parameters areconstrained.
 10. A method as claimed in claim 1 wherein the system ofequations is solved iteratively by minimizing a cost function.
 11. Amethod as claimed in claim 1 wherein the system of equations is solvedusing instrumental variables.
 12. A method as claimed in claim 9 whereinthe system of equations is solved iteratively by minimizing a costfunction.
 13. A method as claimed in claim 9 wherein the number ofequations in the system of equations is reduced by removing equationsfrom the system of equations for which non-zero force terms are present.14. A method as claimed in claim 9 wherein the elasticity of each of theelements is constrained to be between 0 MPa and 1 MPa.
 15. A method asclaimed in claim 9 wherein the viscosity of each of the elements isconstrained to be between 0 kPa·s and 1 kPa·s.
 16. A method as claimedin claim 1 wherein the solving for viscoelastic parameters comprisesgenerating a system of equations comprising a viscosity and elasticityof each of the elements of the model and wherein the system of equationsis treated as being non-linear in the viscosity and elasticity of eachof the elements.
 17. A method as claimed in claim 16 wherein the systemof equations is solved iteratively by minimizing a cost function.18.-24. (canceled)
 25. An apparatus for determining viscoelasticparameters of a tissue, the apparatus comprising: (a) an actuator forgenerating a vibration signal and applying a generated vibration signalto the tissue; (b) a driver circuit connected to drive the actuator togenerate a vibration signal having at least one frequency component; (c)an ultrasound imaging system configured to measure displacements at aplurality of locations within the tissue in response to the vibrationsignal at a plurality of times; and (d) a data processor communicativelycoupled to the ultrasound imaging system and configured to determine theviscoelastic parameters of the tissue including elasticity and viscosityby fitting a finite element model of the tissue to the vibration signaland the measured displacements and solving for viscoelastic parametersof the model, wherein a value for density of each element of the modelis selected and absolute values for one or more of the viscoelasticparameters for each of the elements of the model are solved. 26.(canceled)